/*
* DenseQR.cs
* 
* Copyright (c) 2006, dnAnalytics. All rights reserved.
*/

using System;
using dnAnalytics.Math;
using dnAnalytics.Resources;
using dnAnalytics.Exceptions;

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    /// <summary>
    /// The Compute method locks the object and call InternalCompute. InternalCompute will 
    /// decompose the matrix, overwriting the current matrix with the decomposition, and sets the
    /// the value of q, det and isFullRank.
    /// </summary>
    internal sealed class DenseQR : Algorithm, IQR
    {
        #region Fields
        private bool isFullRank = true;
        private double det = 1;
        private DenseMatrix r;
        private DenseMatrix q;
        private double[] tau;
        #endregion Fields

        #region Constructor
        /// <summary>
        /// Constructs an QR object for the given matrix.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <b>null</b>.</exception>
        public DenseQR(Matrix matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix", Strings.NullParameterException);
            }

            this.r = new DenseMatrix(matrix);
        }

        #endregion Constructor

        #region Public Methods

        /// <summary>
        /// Returns the orthogonal Q matrix.
        /// </summary>
        /// <returns>The orthogonal Q matrix</returns>
        public Matrix Q()
        {
            Compute();
            return q.Clone();
        }

        /// <summary>
        /// Returns the upper triangular factor R.
        /// </summary>
        /// <returns>The upper triangular factor R</returns>
        public Matrix R()
        {
            Compute();
            return r.UpperTriangle();
        }

        /// <summary>
        /// Copies the orthogonal Q matrix into the result matrix.
        /// </summary>
        /// <param name="result">A matrix to copy the orthogonal Q matrix into.</param>
        public void Q(Matrix result)
        {
            Compute();
            q.CopyTo(result);
        }

        /// <summary>
        /// Copies the upper triangular factor R into the result matrix.
        /// </summary>
        /// <param name="result">A matrix to copy upper triangular factor R into.</param>
        public void R(Matrix result)
        {
            Compute();
            r.UpperTriangle(result);
        }

        ///<summary>Calculates the determinant (absolute value) of the matrix.</summary>
        ///<returns>The determinant of the matrix.</returns>
        public double Determinant()
        {
            Compute();
            return det;
        }
        /// <summary>
        /// Determines whether the matrix is full rank or not.
        /// </summary>
        /// <returns><b>true</b> if the matrix is full rank; otherwise <b>false</b>.</returns>        
        public bool IsFullRank()
        {
            Compute();
            return isFullRank;
        }

        /// <summary>
        /// Finds the least squares solution of <b>A*X = B</b>, where <b>m &lt;= n</b>
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
        /// <returns>The left hand side <see cref="Matrix"/>, <b>X</b>.</returns>
        public Matrix Solve(Matrix rightHandSide)
        {
            Matrix result = MatrixBuilder.CreateMatrix(r.Columns, rightHandSide.Columns, MatrixType.Dense);
            Solve(rightHandSide, result);
            return result;
        }

        /// <summary>
        /// Finds the least squares solution of <b>A*X = B</b>, where <b>m &lt;= n</b>
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
        /// <param name="result">The left hand side <see cref="Matrix"/>, <b>X</b>.</param>
        public void Solve(Matrix rightHandSide, Matrix result)
        {
            DenseMatrix rhs = new DenseMatrix(rightHandSide);

            r.lapack.QRSolve(r.Rows, r.Columns, rightHandSide.Columns, q.data, r.data, rhs.data, tau);

            int m = r.Rows;
            int n = r.Columns;
            int bn = rightHandSide.Columns;

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < bn; j++)
                {
                    result.ValueAt(i, j, rhs.data[j * m + i]);
                }
            }
        }

        /// <summary>
        /// Finds the least squares solution of <b>A*x = b</b>, where <b>m &lt;= n</b>
        /// </summary>
        /// <param name="rightHandSide">The right hand side vector, <b>b</b>.</param>
        /// <returns>The left hand side <see cref="Vector"/>, <b>x</b>.</returns>
        public Vector Solve(Vector rightHandSide)
        {
            DenseVector result = (DenseVector)VectorBuilder.CreateVector(r.Columns, VectorType.Dense);
            Solve(rightHandSide, result);
            return result;
        }

        /// <summary>
        /// Finds the least squares solution of <b>A*x = b</b>, where <b>m &lt;= n</b>
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Vector"/>, <b>b</b>.</param>
        /// <param name="result">The left hand side <see cref="Vector"/>, <b>x</b>.</param>
        public void Solve(Vector rightHandSide, Vector result)
        {
            double[] temp = rightHandSide.ToArray();
            r.lapack.QRSolve(r.Rows, r.Columns, 1, q.data, r.data, temp, tau);
            for (int i = 0; i < result.Count; i++)
            {
                result[i] = temp[i];
            }
        }

        #endregion Public Methods

        #region Protected Members
        /// <summary>
        /// Delegates the actual computation to the correct QR subtype.
        /// </summary>
        protected override void InternalCompute()
        {
            q = (DenseMatrix)MatrixBuilder.CreateMatrix(r.Rows, r.Rows, MatrixType.Dense);
            tau = new double[System.Math.Min(r.Rows, r.Columns)];
            r.lapack.QRFactor(r.Rows, r.Columns, q.data, r.data, tau);

            int min = System.Math.Min(r.Rows, r.Columns);

            for (int i = 0; i < min; i++)
            {
                det *= r.ValueAt(i, i);
                if (System.Math.Abs(r.ValueAt(i, i)) <= Constants.AlmostZero)
                {
                    isFullRank = false;
                    det = 0;
                    break;
                }
            }
            det = System.Math.Abs(det);
        }
        #endregion Protected Members
    }
}
